# additional replication code for
# 'Using General Messages to Persuade on a Politicized Scientific Issue'
library(tidyverse)
library(data.table)
library(dtplyr)
library(colorblindr)
library(here)

theme_jg <- function(){
  theme_classic()+
    theme(text = element_text(family = "serif", size = 16),
          strip.text = element_text(face = "bold"),
          plot.title = element_text(face = "bold", size = 20))
}
here::here()
source("code/vaccine_experiments_functions.R")

library(srvyr)
library(grf)
library(multcomp)
library(ggpattern)
library(geomnet)
library(cowplot)

# read in data
wave14.c.exp.oh <- data.table::fread("data/w14_vac_message_cleaned_onehot.csv")

wave14.c.exp.oh$vac_ex_mes_names <- factor(make.names(wave14.c.exp.oh$vac_ex_mes_trunc))
wave14.c.exp.oh$vac_message_resistant <- as.numeric(wave14.c.exp.oh$vac_message == 1)

### grf results

# load up main output
files <- list.files(path = "output/main_output/",
                    pattern = "message*")

for(l in files){
  load(paste0("output/main_output/",l))
}

# join main output
wave14.c.exp.message <- 
  as_tibble(wave14.c.exp.oh) %>%
  left_join(message_harm_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Preventing.harm_binary = predictions_Control_Preventing.harm,
                     variance.estimates_Control_Preventing.harm_binary = variance.estimates_Control_Preventing.harm),
            by = "psid") %>%
  left_join(message_harm_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Preventing.harm_continuous = predictions_Control_Preventing.harm,
                     variance.estimates_Control_Preventing.harm_continuous = variance.estimates_Control_Preventing.harm),
            by = "psid") %>%
  left_join(message_patriotism_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Patriotism_binary = predictions_Control_Patriotism,
                     variance.estimates_Control_Patriotism_binary = variance.estimates_Control_Patriotism),
            by = "psid") %>%
  left_join(message_patriotism_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Patriotism_continuous = predictions_Control_Patriotism,
                     variance.estimates_Control_Patriotism_continuous = variance.estimates_Control_Patriotism),
            by = "psid") %>%
  left_join(message_people_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_People.you.know_binary = predictions_Control_People.you.know,
                     variance.estimates_Control_People.you.know_binary = variance.estimates_Control_People.you.know),
            by = "psid") %>%
  left_join(message_people_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_People.you.know_continuous = predictions_Control_People.you.know,
                     variance.estimates_Control_People.you.know_continuous = variance.estimates_Control_People.you.know),
            by = "psid") %>%
  left_join(message_physician_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Physician.recommend_binary = predictions_Control_Physician.recommend,
                     variance.estimates_Control_Physician.recommend_binary = variance.estimates_Control_Physician.recommend),
            by = "psid") %>%
  left_join(message_physician_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Physician.recommend_continuous = predictions_Control_Physician.recommend,
                     variance.estimates_Control_Physician.recommend_continuous = variance.estimates_Control_Physician.recommend),
            by = "psid") %>%
  left_join(message_scientists_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Scientists.recommend_binary = predictions_Control_Scientists.recommend,
                     variance.estimates_Control_Scientists.recommend_binary = variance.estimates_Control_Scientists.recommend),
            by = "psid") %>%
  left_join(message_scientists_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Scientists.recommend_continuous = predictions_Control_Scientists.recommend,
                     variance.estimates_Control_Scientists.recommend_continuous = variance.estimates_Control_Scientists.recommend),
            by = "psid")

# load up output with facebook covariate
files <- list.files(path = "output/with_facebook_covariate/",
                    pattern = "message*")

for(l in files){
  load(paste0("output/with_facebook_covariate/",l))
}

# join predictions with facebook covariate
wave14.c.exp.message <- 
  as_tibble(wave14.c.exp.message) %>%
  left_join(message_harm_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Preventing.harm_binary = predictions_Control_Preventing.harm,
                     variance.estimates_Control_Preventing.harm_binary = variance.estimates_Control_Preventing.harm),
            by = "psid") %>%
  left_join(message_harm_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Preventing.harm_continuous = predictions_Control_Preventing.harm,
                     variance.estimates_Control_Preventing.harm_continuous = variance.estimates_Control_Preventing.harm),
            by = "psid") %>%
  left_join(message_patriotism_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Patriotism_binary = predictions_Control_Patriotism,
                     variance.estimates_Control_Patriotism_binary = variance.estimates_Control_Patriotism),
            by = "psid") %>%
  left_join(message_patriotism_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Patriotism_continuous = predictions_Control_Patriotism,
                     variance.estimates_Control_Patriotism_continuous = variance.estimates_Control_Patriotism),
            by = "psid") %>%
  left_join(message_people_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_People.you.know_binary = predictions_Control_People.you.know,
                     variance.estimates_Control_People.you.know_binary = variance.estimates_Control_People.you.know),
            by = "psid") %>%
  left_join(message_people_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_People.you.know_continuous = predictions_Control_People.you.know,
                     variance.estimates_Control_People.you.know_continuous = variance.estimates_Control_People.you.know),
            by = "psid") %>%
  left_join(message_physician_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Physician.recommend_binary = predictions_Control_Physician.recommend,
                     variance.estimates_Control_Physician.recommend_binary = variance.estimates_Control_Physician.recommend),
            by = "psid") %>%
  left_join(message_physician_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Physician.recommend_continuous = predictions_Control_Physician.recommend,
                     variance.estimates_Control_Physician.recommend_continuous = variance.estimates_Control_Physician.recommend),
            by = "psid") %>%
  left_join(message_scientists_binary$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Scientists.recommend_binary = predictions_Control_Scientists.recommend,
                     variance.estimates_Control_Scientists.recommend_binary = variance.estimates_Control_Scientists.recommend),
            by = "psid") %>%
  left_join(message_scientists_continuous$predictions %>%
              dplyr::select(1:3) %>%
              rename(predictions_Control_Scientists.recommend_continuous = predictions_Control_Scientists.recommend,
                     variance.estimates_Control_Scientists.recommend_continuous = variance.estimates_Control_Scientists.recommend),
            by = "psid")

# shape
d <- wave14.c.exp.message %>%
  mutate(covnews_fb = ifelse(is.na(covnews_fb), 0, 1)) %>%
  dplyr::select(psid, covnews_fb, starts_with("predictions")) %>%
  reshape2::melt(id.vars = c("psid","covnews_fb")) %>%
  mutate(variable = as.character(variable)) %>%
  mutate(treatment = case_when(
    grepl("Preventing", variable) ~ "Preventing Harm",
    grepl("People", variable) ~ "People You Know",
    grepl("Physician", variable) ~ "Your Physician Recommend",
    grepl("Scientists", variable) ~ "Scientists Recommend",
    grepl("Patriotism", variable) ~ "Patriotism"
  ),
  outcome_type = ifelse(grepl("binary", variable), "binary", "continuous")) %>%
  dplyr::select(-variable) %>%
  group_by(treatment, outcome_type, covnews_fb) %>%
  summarise(gate = mean(value)) %>%
  pivot_wider(names_from = covnews_fb,
              values_from = gate)
names(d) <- c("treatment","outcome_type","no_fb","yes_fb")
d$diff <- with(d, no_fb - yes_fb)

# get prediction comparisons
preds_with_without_fb <- 
  wave14.c.exp.message %>%
  dplyr::select(psid, starts_with("predictions")) %>%
  reshape2::melt(id.vars = c("psid")) %>%
  mutate(variable = as.character(variable)) %>%
  mutate(treatment = case_when(
    grepl("Preventing", variable) ~ "Preventing Harm",
    grepl("People", variable) ~ "People You Know",
    grepl("Physician", variable) ~ "Your Physician Recommend",
    grepl("Scientists", variable) ~ "Scientists Recommend",
    grepl("Patriotism", variable) ~ "Patriotism"
  ),
  outcome_type = ifelse(grepl("binary", variable), "binary", "continuous"),
  facebook = ifelse(substr(variable, start = nchar(variable), stop = nchar(variable)) == "x", "without_facebook", "with_facebook")
  ) %>%
  dplyr::select(-variable) %>%
  pivot_wider(names_from = facebook,
              values_from = value)

# plot predictions
comparisonfb_plot <- 
  preds_with_without_fb %>%
  filter(outcome_type == "binary") %>%
  ggplot(aes(x = without_facebook, y = with_facebook))+
  facet_wrap(~treatment, nrow = 2, ncol = 3)+
  geom_point(alpha = .2, size = .1)+
  guides(alpha = FALSE, size = FALSE)+
  labs(title = "Predicted Effects Comparison, Vaccine Resistance",
       subtitle = "Predicted individual treatment effects with and without inclusion of Facebook use covariate",
       x = "Without Facebook Use",
       y = "With Facebook Use")+
  theme_jg()
ggsave(comparisonfb_plot, file = "results/figures/fb_comparison.png", width = 12, height = 6)

# variable importances

# function to get specific covariate's importance in a specific spec
get_importance <- function(mod, variable = "covnews_fb"){
  df <- mod$importances %>% filter(var == variable)
  names(df) <- c("var","imp")
  return(df)
}

# test
get_importance(mod = message_patriotism_binary, variable = "ideology") %>% 
  mutate(outcome= "Resistance", treatment = "Patriotism")

# write table for importance of covnews_fb across all specs
covnews_imp <- 
  bind_rows(
  get_importance(mod = message_patriotism_binary, variable = "covnews_fb") %>% 
    mutate(outcome= "Resistance", treatment = "Patriotism"),
  get_importance(message_patriotism_continuous, "covnews_fb") %>% 
    mutate(outcome= "Likelihood", treatment = "Patriotism"),
  get_importance(message_harm_binary, "covnews_fb") %>% 
    mutate(outcome= "Resistance", treatment = "Preventing Harm"),
  get_importance(message_harm_continuous, "covnews_fb") %>% 
    mutate(outcome= "Likelihood", treatment = "Preventing Harm"),
  get_importance(message_scientists_binary, "covnews_fb") %>% 
    mutate(outcome= "Resistance", treatment = "Scientists Recommend"),
  get_importance(message_scientists_continuous, "covnews_fb") %>% 
    mutate(outcome= "Likelihood", treatment = "Scientists Recommend"),
  get_importance(message_physician_binary, "covnews_fb") %>% 
    mutate(outcome= "Resistance", treatment = "Your Physician Recommend"),
  get_importance(message_physician_continuous, "covnews_fb") %>% 
    mutate(outcome= "Likelihood", treatment = "Your Physician Recommend"),
  get_importance(message_people_binary, "covnews_fb") %>% 
    mutate(outcome= "Resistance", treatment = "People You Know"),
  get_importance(message_people_continuous, "covnews_fb") %>% 
    mutate(outcome= "Likelihood", treatment = "People You Know")
) %>% dplyr::select(treatment, outcome, imp) 

covnews_imp %>%
  write.csv(file = "results/tables/covnews_imp.csv")

covnews_imp %>%
  knitr::kable()
